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NOMENCL^^TURE 


Greel^ Symbols: 

= 1 for Fiber 

^ =2 for Matrix 

= Small time interval 

= Small temperature difference 

fll / n 2 n -5 = Natural coordinates for the triangular 

finite element 

© = Angular coordinate 

0 = Field variable 

4 = Specific heat per unit volume 


Alphabetical Symbols: 


A 



D 

f 


[P] 

g 

h 

K 

[k] 


= Area of the 2-D finite element 
= Constants 
= Constants 

= Domain of the finite element 
= See table 2.1 
= Heat flux matrix 
ss See table 2.1 
= See table 2.1 
= Thermal conductivity 
= See table 2.1 
= Thermal stiffness matrix 
= Interpolating function 
= Radial coordinate 



0 


"0 


/ 


X 


= Radius of fiber 
= Radius of equivalent cell 
= Domain boundary segments 
= Temperature 
= Applied temperature 
= Time 

= Impulse duration 
= Axial coordinate 


Superscripts: 

a , (a) = a - constituent 

(e) = Finite element (e) 


Subscripts : 


1 

2 

X 

r 


= Fiber 
= Matrix 

= In the axial direction 
=: In the radial direction 


i # j / h 


= Nodes of triangular finite element. 



ABSTRACT 


A formulation has been developed to study the 
phenomenon of one-dimensional Heat Conduction through 
Unidirectional, Continuous Fiber Reinforced Composites. 
Further, a transient Finite Element procedure has been 
developed ho solve the above formulation for variations 
in thermal properties and composite geometries - 

An initial boundary value problem has been 
selected and solved for temperature-dependent thermal 
properties. System sensitivity to non-linearities aris- 
ing from temperature dependence of thermal properties 
has been examined extensively for composites having a 
constant Piber»-Matrix volume fraction. 

The Finite Element Method results have been 
presented in the graphical form and discussed. It has 
been found and demonstrated that non-linear behaviour of 
thermal properties could significantly alter system 


responses 
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Chapter 1 


INTRODUCTION 


1.1 COMPOSITE MATERIALS 

A composite material, in general, could be defined, 
as one having two or more distinct constituent materials 
or phases. They offer outstanding mechanical properties 
coupled with flexibility in design and ease of fabrica- 
tion. The Use of composite materials has become quite 
common, ' these days, in industries like aircraft, auto- 
mobile, sports-goods and electronics. Notable among the 
many composite materials, are the fiber composites which 
have the ability to make use of the outstanding strength, 
stiffness and low specific gravity of fibers such as glass, 
graphite or Kevlar. Also, in the case of fiber composites, 
the ability to synthesize the material structure during 
manufacture is a great advantage. 

1.1.1 Characteristics 

Composites consist of one or more discontinuous 
phases embedded in a continuous phase called the matrix. 

The discontinuous phase is usually harder, stronger and 
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primarily meant for reinforcement. An exception is the 
class of rubber-modified polymers, which have a rigid 
polymer matrix filled with rubber particles. 

The properties of composite materials are strongly 
influenced by the properties of their constituents, their 
distribution and the interaction between them. They also 
depend on the geometry of the reinforcement with reference 
to the system. Concentration , measured in terms of volume 
or weight fraction, is perhaps the parameter with greatest 
influence on composite properties. It being controllable 
during manufacture, desired properties of the composite 
could be created with great facility. The concentration 
distribution of the particles, which depends on their 
spatial relations, is a measure of the homogenity of the 
system. The orientation of the reinforcement affects the 
isotropy of the system. 

7.1.2 Classif ication 

The strengthening mechanism of a composite, 
strongly, depends on the geometry of the reinforcement 
and hence, it would be useful to classify them accordingly. 

Composite materials can be, in this light, 
broadly classified as fibrous composites and particulate 
composites. Fig. 1.1 gives a more detailed classifi- 
cation. 



COMPOSITE MATERIALS 


FIBER REINFORCED PARTICLE REINFORCED 

(FIBROUS COMPOSITES) (PARTICULATE COMPOSITE) 
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] 
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The chief characteristic of a particle is that 
it is non-fibrous in that, its length is not very much 
greater than its cross-sectional dimensions. As our study 
is related to fibrous composites, let us have a more 
detailed look at their characteristics, 

1 4 1 . 3 Fibrous Composites 

With normal engineering materials, imperfections 
in a direction normal to applied loads are particularly 
detrimental to their behaviour. Man-made fibres have 
considerably lesser flaws perpendicular to their lengths 
and are therefore, stronger, resulting in their use as 
reinforcements . 

The most important reinforcement fiber is E-glass 
because of its relatively lower cost. Boron, graphite 
and aramid polymer fibers (Kevlar 49) are useful for their 
high stiffness. 

Fibers, because of their small cross-sectional 
dimensions, are not directly utilized in engineering appli- 
cations but are used as fibrous composites by embedding 
them in matrix materials. The matrix serves to bind the 
fibers together, transfer loads to the fibers, and pro- 
tect them against environmental attack and damage due 


to handling 
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Fibrous composites can be broadly classified as 
single-layer and multi-layer (angle-ply) composites. 

Single-laver " composites may/ actually, be made-up of 
several distinct layers, each layer having the same orien- 
tation and properties. The entire laminate, thus, behaves 
like a " Single-layer" composite. "Multi-layer" compo- 
sites, on the other hand, consists of several layers, 
each layer or lamina being a single-layer composite. The 
orientation of these single-layers is varied according 
to the desired design. 

When the constituent materials in each layer are 
the same, they are called laminates. Hybrid laminates 
are multi-layered composites with layers made up of diffe- 
rent materials. 

The chief advantages of employing fibrous compo- 
sites lie in their high strength-to-weight ratio and their 
controlled anisotropy which facilitates variation in 
desired ratios of properties in different directions. 

1.1.4 Fibrous Composites as Insulators 

The combination of high stiffness and strength 
with low weight enables fibrous composites to be natural 
candidates for aerospace applications. They are employed 
in supersonic transport to accommodate thermal stresses. 
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in sppce transport as heat shields and nose~tip insula- 
tors, and in space structures to reduce thermal deflec- 
tions . 

An outstanding example of their use as insulators 
is the case of the reusable surface insulation (rsi) tiles 
employed on the space shuttle orbiter [l]. Their cost 
effectiveness, resistance to thermal shoch and high insu- 
lative efficiency mahe them the best choice for such 
applications. The RSI tiles consist of silicon fibers 
dispersed in silicon carbide. 

Heat treated polyacrylonitrile (PAN) based carbon 
felt, in a rough laminar carbon matrix, gives a material 
with high resistance to thermal shock [ 2 ]. PTiN-felt carbon, 
in chemically vapor deposited carbon matrices, also are 
good insulators [ 3 ]. 

1.1.5 Heat Diffusion in Piber Reinforced Composites 

It has been observed that temperature and mois- 
ture are the primary environmental conditions that affect 
the behaviour of a structural composite during its service 
life. Hence, the analysis of thermal and/or moisture 
fields in a composite is of great importance. 

Conduction, parallel to the fibers of an uni- 
directional fibrous composite, is a problem of conside- 
rable practical importance for those cases in which such 
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composites are used as thermal insulators. Re-entry 
vehicle heat shields and nose-tios, in which the primary 
direction of heat flow is the fiber axis, are typical 
examples . 

The reduction of these problems to the one dimen- 
sional type entails greater computational efficiency. Most 
materials, used currently for these purposes, exhibit 
temperature dependent diffusive properties. The result- 
ing nonlinearities render analytical treatment intract- 
able. lA/hile considerable work has been carried out on 
the finite element analysis of conduction, the studies 
were limited to linear cases [ 4 ]. It is, therefore, of 
interest to analyse the non-linear heat transfer under 
such circumstances. 

ii2 FINITE ELEMENT METHOD IN THERMAL ANALYSIS 

The Finite Element Method (FEM) was first applied 
to structural analysis. The approach is now recognised 
as applicable to any problem that can be formulated in 
the differential or variational form. This recognition 
is a boon to the analysis of thermal fields and coupled 
problems which arise in many applications. Aerospace 
structures, specially, can operate in severe thermal 
environs which significantly affect their structural 
design. Examples include convectively cooled structures 
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under development for hypersonic transport, thermal pro- 
tection systems for the space shuttle and large flexible 
space structures currently being designed for earth orbit. 
The design of these systems requires a strong interaction 
between thermal and structural analyses. The FEM offers 
the greatest potential for efficient coupling of the two 

[5]. 


1.2,1 Essence of the Method 

The heat diffusion equation could be written as, 


v^kv't + q_c|| = o (1.1) 

which is a classical example of a continuum mathematical 
problem. 


The numerical 'discretization' of 
nuum problem with its boundary conditions 

T - T = 0 on 

n'^ 1< V T _ q 0 on 


such a conti- 


( 1 . 2 ) 


is the essence of all numerical solutions, where and 
9 ue portions of the boundary on which temperatures T 
or fluxes q are given and n represents the unit normal 
vector to the boundary. 

The finite element method can, in general, be 
defined as a process wherein. 



9 


(i) the continuous function T is approximated by 
a series of parameters and specified trial 
functions N^(x,y, 3 ) in the problem domain Q as 

T,^:?T = 2N^a^; i = l...n (1.3) 

(ii) writing the differential equation (such as (1.1)) 
and its boundary conditions (such as (1.2)) in 
the general form; 

A(T) = 0 in Q 

(1.4) 

B(T) =0 in w 


The approximating equations, from which the 
solution is to be obtained, are written as a set 

/_ W. A ( T ) dQ + / W. B ( T ) dm = 0 
£3 J to 1 

J = 1 n (1,5) 

where Wj and Wj are suitable weighting functions which 
ensure that as n 


A ( T ) -» 

and B ( T ) — <► 

i.e., at all points, 
the exact solution. 


0 in Q 


( 1 . 6 ) 

0 in CO 

the approximating solution tends to 


The definite integral of equation (1.5) can be 
written as simply a sum of such integrals occurring on 
'subdomains' into which the whole 'domain' is divided 
(Fig. 1.2). Thus, if 




. :nm:is2 ami%aocHMg ^ maasm 
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p. a 

Q = S Q and w = S w (1.7) 

0=rl 

then, for all finite functions ( ) / we have, 

m 

/^ ( ) dQ = S ; ^ ( ) dQ 
S3 . — e 

e=l Q (l^g) 

m 

/ ()da)= S 

e=l w 

where, and are associated with 'elements' into 

which we divide the problem. 


This allows the whole region to be divided into 
standard type of sub-regions, in which the parameters 
aj^ are usually the nodal values of the independent func- 
tion T and in which the tr^al functions are defined in a 
local manner. The integrals can, then, be evaluated 
element by element and the approximating equations (such 
as (1.5)) obtained by a simple addition of element con- 


tributions . 
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1.3 PREVIEW OF PAST WORK 

On reviewing the work done in the area of con- 
duction, it was found that while literature is rich with 
various methods for the solution of linear problems, there 
is relatively little material available on non-linear 
problems . 

In 1951, Storm [s] showed that the mathematical 
condition for conversion of the 1-D non-linear conduction 
problem to the linear case is the constancy of [l/(KS)^'^^] 
[log (S/K)^^^]. The conversion, however, was found 
valid only for simple metals. 

Pattle [ 7 ], in 1959, presented a closed form 
solution to the non-linear problem using an integral 
method for the restrictive case of an internal heat source 
in an infinite medium whose thermal conductivity follows 
some power law K(t)''^t’^. During the same year, Yang et al., 
[s], discussed the integral method for multi con^onent 
planar solids, 

Tsang [ 9 ], Andre-Talamon [lo] and Olsson [ll], 
in later years, discussed small perturbation methods to 
solve non-linear conduction problems. However, these 
methods are complicated to use, as it becomes difficult 
to obtain solutions to higher order perturbation terns. 
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Goodman [ 12 ] and Imber [is] have described inte- 
gral or heat balance methods of considerably low accuracy. 
The latter amended the integral method to include phase 
change problems. 

Padovan [l4l/ in 1974 used complex series repre- 
sentations and properties of adjoint differential opera- 
tors to develop a finite element procedure for non-linear 
conduction, which is rather complicated to use. 

Leelamma et al./ [ 4 ], in 19S1, have pointed out 
the non-availability of finite element solutions to non- 
linear conduction problems, thus far, and proposed a per- 
turbation method for weakly non-linear cases of an isotro- 
pic material. 

Analysis of behaviour of composite materials 
began, rather, late. It was only in the early seventies 
( ^Q 7 ;[97 3 ) that some research was carried out on diffu- 

sion in composites. Bedford and Stem [l5,l6] proposed 
a continuum theory for diffusion in composites, modelling 
the constituents as superimposed continue. The same 
theory was found effective for one-dimensional wave pro- 
pagation in composites in the direction of reinforcement. 
Hegemier and others [l7,13] presented a continuum theory 
and an analvtical procedure for modelling wave propaga- 
tion in laminated and uni-directional fiber composites. 
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The theory takes care of mixture interaction and micro- 
structural information, but is restricted in application 
to wave propagation only. 

Horvay et al#, [l9]/ first studied the case of 
transient heat conduction in laminated composites , 
omitting the case of fiber composites. Padovan [soj 
gave an analytical solution for transient temperature 
distribution of a 2—D general anisotropic half space. 
Later, using a variational procedure, he developed a 
finite element approximation for steady heat conduotion 
in anisotropic media [21,2?]. The methods, however, are 
not applicable for transient and non-linear situations. 

Haltman et al., [ 23 ], examined experimentally 
and numerically, steady heat conduct ion in a two phase , 
two-dimensional model of aluminium bars dispersed in 
foamed polyurethrane. Hegemier et al., [24,25,26], gave 
various continuum theories related to fiber reinforced 
composites and laminated composites but, applicable to 
the phenomenon of wave propagation only. In 1975 ,lSIayf eh 
[ 27 ] gave a continuum mixture theory for heat conduction 
in laminated wave guides, using approximations enabling 
analytical solutions . 

Chen et al . , [2s], calculated the thermal con- 
ductivity of composite materials containing thinly 
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dispersed spheroidal inclusions. The associated problems 
of heat flow across one such inclusion embedded in an 
infinite matrix were solved, in the presence of thermal 
boundary resistance. Maewal et al . , [29], developed a 
continuum model for wave propagation problems, and applied 
it for diffusion type processes also, in a laminated 
composite with periodic micro structure and obtained 
finite difference solutions. Again, Maewal et al., [30] 
and Murakami et al., [31] developed mixture theories for 
the process of heat transfer in unidirectional fibrous 
composites with periodic hexagonal microstructurek 

Patankar [ 32 ], gave a finite difference solution 
of heat conduction in the case of non-uniform thermal 
conductivity in a composite slab. Poon [33], reviews 
the various methods to solve problems of heat conduction 
in laminated composites. Murakami et al., [34], extended 
their earlier work [30], to include non-linear material 
behaviour. Laura [ss], gave analytical and numerical 
solutions to determine the temperature distribution in a 
composite rod with variable internal heat generation. 

It was in April 1983 , that Arun Aggarwal [ 36 ] 
developed a Finite Element Analysis package for the 
longitudinal heat diffusion in fiber-reinforced compo- 
sites. The work, however did not consider non-lineari- 
ties arising due to temperature dependent material 
,properties.^ : 
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The Present Work 

In the present work, non-linearities arising due 
to temperature dependence of material properties have 
been taken care of as an extension of the work reported 
in [ 35 ]. An attempt has been made to show how non- 
linearities could affect the heat diffusion and tempera- 
ture distributions. For this purpose, the author deve- 
loped a Finite Element software package and executed it 
on the DEC 1090 system. 

An initial boundary value problem was selected 
and solved for non-linear material behaviour. Polynomials 
were used to represent the nature of non-linearity and 
the entire excercise was carried out for two different 
sets of materials. It was found that non-linearities 


could significantly alter the temperature- distributions 



Chapter 2 


FORMULATION AND ANALYSIS 


2 . 1 THE FIELD PROBLEM 

A large class of field problems are governed by 
similar partial differential equations for the field 
variable 0. Heat conduction, wave propagation and con- 
vective ditfusion are three examoles. 

Consider that it is desired to evaluate the 
field variable 0 in a three dimensional solution domain 
Q, bounded by surface S (Fig. 2.1). The field equation 
to be solved could, in general, be written as, 


d0 


ftx + S? SI- SS 


^0 


f(x,Y,z) 

( 2 . 1 ) 

where K , K , K and f are given functions. The phy- 
X y z 

sical interpretation of the parameters in the above equa- 
tion depends on the particular physical problem. Table 
(2.1) lists some typical problems and their parameters. 

The boundary conditions of (2.1) could be written 

as , 


0 


0(x,y/2;) on S 


1 


(2 . 2 ) 




Fig. 2«1 - 3-D SOLUTION DOMAN 



No. Problem Field , K , 

parameter ^ 

0 ^2 


f g b 


1 . 


Dif fusion~f low Hydraulic 
in porous head 

media 


Hydraulic Internal Boundary - 
conductivity source flow 
flow 


2. Electric 

conduction 


Voltage 


Electric Internal Exter- 

conductivity current nally 

source applied 
boundary 
current 


3. Electro-static Electric Permitivity 
field force 

field 

intensity 


Internal 

current 

source 


4. Heat 

conduction 


Tempera- 

ture 


Thermal Internal Boundary 

conductivity heat heat 

genera- genera- 
tion tion 


Convec- 

tive 

heat 

trans- 

fer 

coeff i- 




5. Torsion 


6. Seepage 


Stress Reciprocal 

function of shear 
modulus 


Angle of 
twist 
per unit 
length 


Pressure 


Permea- 

bility 


Internal 

flow 

source 


7. Magnetostatics Magno- 

motive 

force 


Magnetic 

permea- 

bility 


Internal Exter- 
magnetic nally 
field applied 
source magne- 
tic 
field 
inten- 
sity 


TABLE 2.1 : I DENTIFICATION OF PHYSICAL PARAMETERS 



20 


and 


K 


X 




n 

X 


+ K 


n 

Y by Y 


+ K 


^0 


2 d2 


n 

z 


+ g (x/Y/z) + h (x/V/z) = 0 on S 2 (2.3) 

where g and h are known apriori and n^ , n^ , are the 
direction cosines of the outwar^^ normal to the surface. 


2 .2 FORMULATION 
2.2.1 Composition 

Consider a two phase material whose phases are 
of cylindrical shape, with all generators oriented in 
the same direction. For convenience, let us concern 
ourselves with materials of cylindrical shape with gene- 
rators parallel to the phase region generators. It is 
also assumed that phase cylindrical regions are unin- 
terrupted along their lengths. 

In the general fiber reinforced material, there 
is no specific geometrical distinction between the two 
phases. If we consider one of the phases as consisting 
of cylinders embedded in another phase, we will call the 
embedded cylinders the FIBRES and the phase in which they 
are embedded, the MATRIX. 
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2.2.2 Georoe'bry 

If we assume the fibers to be of a specific shape 
(in this case circular, say), the geometry is one of regu- 
lar arrays of identical fibers in a matrix. The most 
general binds of arrays are perhaps the square, rectan- 
gular and hexagonal arrays. 

For the present study, we consider a hexagonal 
array of circular, cylindrical fibers ( constituent- 1) in 
a matrix (constituent-2), as shown in Fig. 2.2. 

Using cylindrical polar coordinates, the domain 
could be described as occupying a space 0 ^ ' 

(Fig. 2,3). 

2.2.3 Foimulation 

We consider the temperature or heat flux boundary 
conditions at x = 0 and at x = L to be such that the tem- 
perature distribution is similar in each hexagonal cell. 

As a consequence, the external boundary of the unit cell 
becomes a line of symmetry, and the component of the heat 
flux normal to the boundary vanishes. It is further 
assumed that there is no thermal resistance across the 
fiber-matrix Interface. 

To describe the heat conduction in a typical 
unit cell, the hexagonal boundary is approximated by a 






Fig. 22 - HEXAGONAL ARRAY 


COMroSITE CYLNtXR CONSTRUCTION 






I 


li 


Fig. 2-3 ' GECMiTRY AND COCHDINATE SYSTEM 



circle 


2 4' 


. As a result, the temperature distribution within 
the cell is axisymmetric with zero heat flux normal to 
the outer circular boundary of the matrix. 


The initial boundary value problem that des- 
cribes transient heat conduction in the cell is then 
given by the following equations: 


(a) Conservation of energy 


1 

r 


d 

dr 


(r q. 


(r^) 


A- 

dx 




(a) 


4 


(a) ^ 

dt 


(a) 





(2.4) 

(b) 

Constitutive equations 



q <"■> 

_ _ K 

X 6x 

(2.5a) 


(„.) 

^r 

(a) 

== ~ r 2>r 

(2 .5b) 

( c) 

Interface 

conditions 



jCl.) 

(x,r^ , t) = 

(2.6a) 



(x, r^ , t) = 

( 2 . 6b ) 

(d) 

Symmetry 

condition 




CO 

ft 

il 

O 

(2.7) 

(e) 

Initial and boundary data. 



. . (a) (a) 

The quantities q^ , q^ ' ^ 
and K in the above equations are, respectively. 


the heat flux components, heat capacity, temperature and 
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thermal conc^uctivity components of the oc ( = 1 or 2) 
constituent. 

The problem has been solved to determine the 
evolution of the temperature field in the half space x > 0 
subject to a step function at the boundary/ 

T (O, r, t) = T^ , 0 < t < t^ 

(2.3) 

0 , t > t^ 


2.3 FINITE ELEMENT FORMULATION 


The governing equations for the axi symmetric 
heat conduction problem could be written as 


1 A (K r i?) f 


(K 


^T 


) 


4 




(2.9) 


b dr ^ dx ^^x ax' - -- at 

along with its boundary conditions (Fig. 2.4). 

As we consider the temperature evolution in the 
half space x > 0 , the boundary conditions at the two ends 

become 


At X = 0 


At X = L 


T = T 


0 ' 


0 <_ t <_ tg 


( 2 . 10 ) 


id 

dx 


0 


= 0 


t > t, 


'0 


By symmetry of the composite cylinder assembl- 
age, the heat flux normal to the cylindrical surface of 

the unit cell vanishes. Thus, 

'dT _ Q (2.11) 

at r- = r 2 ^ - 0 
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Due to the symrnetry of the geometry and boundary 
conditions, it would suffice, if analysis were carried out 
for one half of the cylindrical cross-section. Thus, due 
to symmetry , 

at r = 0 , = 0 and T is finite (2.12) 

dr 

We shall resort to the Galerkin's approach to 
derive the FEM equations. 


2.3.1 Finite Element Discretization 

without deciding about the type of element to 
use, we can write 

(x, r, t) = U,rf {TCt)! 

(2.13) 

where, the interpolating functions need to preserve 
C° continuity only i.e., continuity of the field variable 
at the element interfaces. 


The Galerkin's criterion could then be applied 


to obtain 


^(e) h [ ^ 


r 


) + ^ (K 


Sx X dx 


St 


4 ^ ] dD^'"' = O 


.(e) 


(2,14) 


,(e) . 


where i = l,2,...n and D^ is the domain of element e . 


For axisymmetry/ 

dD*'®^ = 2 "^r dA^®^ 


(2.15) 



28 


Using this in (2.14), we get, 


_ N . [ i (K_r 1-?-^) + ^ (K, 


(e) 


D 


( e) 


i ^ r dr ^ dr 




( e) 


- 4 


dx X ax 

11^] 2 - r aA‘-> = 0 


(2 . 16) 


Integrating, by parts, the first term w.r.t. 'r* and the 
second term w.r.t. 'x', combining and rearranging the 
terms, one obtains 


dT 


( e) 


;; [n, ^ (K„r ^ ) + N, „ (K.r 

( e) 


i/r r dr 


D 


i,x X ax 

XT' ^ ® ^ ( p,) 

+ ;/ 4 r N. — 


a 4 e) 


D 


(e) 


i dt 


X 


J 


top 


N. (K r 


%^(e) I ^right 


bottom 


1 r dr 


dx + 


'left 


^right 

+ J N, (K .r gi ) 

^lef t 


■ \ AX. .. 

1 X dx 


X 


top 


X 


dr 


bottom 


or. 


(2.17) 


uCe) 


i,x X ax 
,(e) 


(e) 


jzT (K^r 

.(e) 


i 5^ ’"r + "^x^ dSE ^x^^^ 


(2 .18) 
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Introducing the boundary condition for flux on the element 
( e) 

boundaries ^2 , we havr 

K r ^ n + K r ^ n + r q (t) = 0 for t > O 

r r X Qx x ^ 


(2 . 19 ) 


Substituting in (2. IS), we get 


SI [ + 


Ae) 


+ ^ r Ni 

T^(e) 


dA*-''^ + jZT' N. r q dS^ = 0 

<,(e) ^ 

^2 

( 2 . 20 ) 


Using from (2.13) in (2.20), 

;; [K^r 1 N/r;^ + K^r :N,x|'^] dA^®^ + 


+ JS [nr N. dA^®^] [t] 

(a') 


^ jf N. r q dS 2 ^^= 0 


(2.21) 


where T = 


the time derivative of the temperature. 


Casting (2.21) into the matrix form, we obtain. 


S! (K r 

-r~>( ® ) 


{s,r] |N,r}^ + V !n,x| 


+ // 


(U r 


nU csa'®’) qd"' 


+ ^ iN^ r q ds)' 


= 0 


( 2 . 22 ) 
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+ [P 


This is the basic element matrix equation 

For each element, 

the thermal stiffness matrix. 


; ( el 


= JJ (k r ^N,r ] /ibri^ 


D 


<e) J 


+ 




{n,x 1, [M,xr ) 

‘ J I I 

the thermal capacitance matrix 




[c] 


( e ) 


= JS 4 r I N i 'n|'^ dA^®^ 
(e) ’ J I ^ 


D 


and the heat flux matrix 
[p]^®^ = jzT 


(e) 


j N "j* r q*'® ^ 


0 (2.23) 


(2.24) 


(2.25) 


(2,26) 


Our problem is an axisymmetric case, which 
requires continuity of the field variable alone. Let us 
choose triangular ring elements (Figs. 2.5, 2.6) with 
linear variation of field variable T along the element 
boundaries as it satisfies both C° compatibility and 
completeness requirements. For the three node element 
we have. 


T 


(e) 


T 

pi) ; 

' T^(t) ' 

; n2 V 

TjCt) 1 

1 Hs! 

^TjCt) j 


(x,r, t) 


(2.27) 
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For linear variation of T along the element boundaries ^ 

„ H 

1 1 “ ngjj - rjX.)+(r.x^ - r^x'j+lrj^x,- r.xj 

The denominator of is equal to 2 a where A is 
the area of the triangleo Likewise/ it is also the deno- 
minator of ri 2 '13* 


H ence , 


'^2 ~ ^2 


(rix - r X. )-)-(r x^ - r.^x) + (r^Xi - r^^x^) 


^ 3 ” ^3 


- r.x. ) + (r . X - r x. - r. x ) 


(2.28) 


Also / 

+ TI2 + TI3 = + N2 + = 1 

We have as a consequence of the above definitions. 



^ - 


fl’ = 


1 ,r 

“ ar 

2A 


2A 



(>'1; - 

li’ = 

®2 

2 ,r 

a r 

2A 


2A 



(Xi - 


®3 

3,r 

ar 

2A 


2A 


where , B2 and are constants 


(2.29) 




*'’1 _ - ■ 

4,X “ BX “ 2A 

‘=^1 - fa 

^9 .-sf " 2 A~ 2A 


L ihewise 
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(r . - r. ) c. 
N = — J-- ^ _ _3 

3fX 2A 2A 


where , C 2 and are constants. 

We now proceed to evaluate the elemental 

matrices . 

'thermal conductivity (s tiffness) matrix: 

We have 

( p ) 


(2.30) 


[k] = // ^N,r^ / N/r^-"^ 

jj(e) ( J ^ 


+ K.r jN,x| fN,x\'^)dA 


( e) 


(2.24) 


X- L ■- j [ ' -j 

Our thermal conductivity K is temperature depen- 
dent for its value. As suggested by Abel and Desai [46], 
assuming piecewise linearity for K i.e., a constant but 
different thermal conductivity for each increment, we 
could proceed as though the case were linear. Using this 
ideology and substituting for N, r and N,x from (2.29) 
and ( 2 , 30 ) we obta in , 

T 


[k] 


k r 
r 





f B, / 



1 ; 

®2 i ^ 

' 



B 1 

1 3 J 

3 j 


K r 

X 


4A 


(e)- 



( 


1 


( 

1 C 

‘1 


V 3 

] 

1 ^ ) 


dA 


( e) 


(2.31) 
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AIq 6*13178 1 C triRnlpulB t ion will pirovB ir = ir.r] 4-17*1^ + 3 rT) 

1 1 3 2 Ic 

Hence, 


[k]<' 


I< 


// , 


BJ+C^ BjB^+CjCj BjB 3 + CjC^ 


Symm , 


2 2 
B^+C 

2 


^2^^+ C2C2 

2 2 i 

J 


(Ti + r, ri2 + 113) dA 


( e) 


.... (2.32) 

For a triangular element, the Integration formula is 


a B r (p.) ) 

/ ri%nhodA^^^= 2A 

^(e) ^ ^ (a+P+r+2)i 


(2.33) 


Employing this, (2.32) transforms to 

>? 


[k] 


( e ) K 


(e) 


4A 


(e)' 


shc^ BjBj+CjCj B^B3 +CjC3' 


9 2 

^2+^2 


Symm. 


®2®3‘^^2^3 


J 


tJ (r^Bi+r'n2+>^V:n3)aA 


(e) 


D 


[i<] 


(e) K 


(e) 


(e) 

bI+cI B^Bj+CjCj B3B3+C3C3- 


a2A 


(e)' 


2 - 2 


(ri+rj+r^^) 




2 2 
^3-^S 


(2.34) 
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, ((=>■) 

In equations (2.26) and (2.28)/ K has been 
used ns IF it were a constant and eaual in both axial and 
radial directions. This has been done for convenience of 
expression. During solution, it is necessary to club the 
value of radial conductivity with the 'B" terms and axial 
conductivity with the 'C' terms. 


( i i ) Thermal capacitance matrix ; 
We have 


[c] 


r Jn i I d-A^®^ 

(e) Ml] 


( 2 # 2 ) 


Extending the argument of piecewise linearity to the 
thermal capacitance and considering it constant over an 
element, we have 




(e) 


D 


// 


;T 


+ n3)dA 


(e) 


ti 2 n n^, h, 

112 1 , 


'2 2 3 


h. 


(e) 


Ui q + g h + 

Employing the Integration formula, this refluces 
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^ 60 

( 6r. +2r .+2r. ) (2r. +2r.+r, ) (2r. +r.+2r, ) 

IJK ijk IJK 

( 2r. +6r .+2rn, ) (r. +2r .+2r. ) 

IJK IJK 

Symm. {2r^+2r ^+6r^) 

. o . , « (2 . 35 ) 


( i ii ) Element heat flux vector ; 
From (2.26), 



This term exists only for elements at the boun- 
dary where the heat flux is specified. Oonsider that on 
the boundary edge the heat flux is constant for any ele- 
ment, This assumption simplifies the analysis. 

Thus , 

' ^ V e / 

^2 



(t) 





Ui *^1 *^2 J 



For a linear element, the appropriate integration formula 
is 
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P! 

("+P+1 }! 


whcro is the length of the linear element. 

Using this, (?.36) transforms to. 



S, 


r 

-n 


r ./3 4- r ./6 


r./6 + rVB 


(2. 37) 


where S^__j is the length of the boundary edge of the 
element. 


2.3.2 General Matrix E qu ations 

The element matrices as given above are calculated 
for each element and assembled together to give the system 
of simultaneous linear equations in the matrix form as 

[k] I t| + [c] ( T I + [f} = 0 (2.38) 

Our problem being a transient case, we are con- 
fronted with time derivatives of temperatures. Zienkiwicz 
[ 37 ] shows that either finite difference aonroximations 
or finite elements in time (use a shape function defined 
in the time domain) may be used for computer solution 
and that the methods are equivalent. 

We will adopt the former method, replacing time 
derivatives by finite difference approximations. 

Prom (3.38) at time t, 
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[^t] + [c. ] 




l"t] 


( 2 . 39 ) 


At time t“5t 


[’'t-et^ bt-6tl + f'=t-6t] ^ ['’t-6tj = ° 


( 2 . 40 ) 


Then, 


+ it ) T ( , 6_t m / 

J t-6tj + 2 I ^tj + 2 ( \-6tj 


f^t] = ^ ■ Tt-st} >- iit-et] 

Substituting equation (2.41) in (2.?9) we have, 
^’^t^ ^'^t^ ‘'5t^ I'^t] ~ ["^t-St^ ^ 


1 - = 0 


([k^] + ^ 


F rom ( 2 , 40 ) , 


[c^]) = [cj t ; Vat} + It K-St} 


r 

- 


(2.42) 


It, = [c^ (-[k, ,.] It,. .[ - If, ) 

I, t“6t] ^ t-et"^ ^ t-6t ]_ t-6tj (_ t-6tj 


(2.4?) 


Substituting this in (2.42), we obtain. 


([k.] + I 


t-' ^ [Ct])[T^| = [c^](-[c^_g^] 


Rearranging, 
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([i<p + It igh = ( Ittcp- 


- fq] ht-6tt 

- [Ct] 

(2.44) 


For most materials the variation of thermal conductivity 
with temperature is much more pronounced than that of 
the thermal capacity. Considering, therefore, that this 
variation be neglected in comnarison, we can write 

[c^_ J = [C,] = [c] (2.«) 

-•-([k^] + It [C]) = (ft [C] - [Kt.ith 

- ( pti 

This equation is valid for transient analysis. 

For steady state problems, the time derivative of tempe- 
rature is zero, reducing (2.38) to the form of the Fourier 
law of conduction, 

[k] T + I F^ = 0 (2.47) 

solution oF this equation gives the steady state results 
in one step of computation. 

2.4 SOLUTION PACKAGE 

A computer program has been formulated to solve 
the matrix equations in FORTRAN. The program enables 
holding Of boundary nodes at constant temperatures. 
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As we arc solving for the temperature field, it would 
be of interest to suppress the nodes for which the tem- 
perature is already known. 

We have the eguation 

[C] - [K^. J] 

- ‘ ihl + K-60 > (2.46) 

The solution procedure, essentially iterative, 
is as follows; 

The initial temperature field being specified, 
the value of conductivity matrix elements is found using 
the specified conductivity temperature relations . U sing 
those values of conductivity , the equation (2.46) is 
solved to obtain the temperature field at the end of the 
time interval fit. The values accruing from the new 
temperature field are then the basis for calculation of 
conductivity values to be used in the computation of the 
temperature field during the succeeding iteration. The 
process is repeated for desired number of iterations or 
for tho desired time period. 

The matrix equation (2.46) can be rewritten to 
enable suppression of nodal equations with constant or 
fixed temperatures. Denoting the conductivity matrix 
by [a] and the heat capacity matrix by [b] we can 


write. 
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[a+bJ^tI = [b-a] |T^| (2.48) 

where, solution for the vector T is sought and vector 
■ Tq is the temperature field, just computed or at the 
beginning of the iberation. If ' T represents the fixed 
temperatures, ('2.43) can be rewritten as. 





which leads to. 



C-hi + hiMh) = - '"o) 

- 2 Ka + ha] fhi] [h} 

+ ? [b^j] [ Tj 

= - atijj] - [Ajj + 

• ' • ]-*ll ®ll] {’’ll = ' 2 [*12-] 1^1 ]-®ll"*ll] (hj 

(2.49) 

In the solution of the above eguation it would 
be of considerably economy if use were made of the 
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symmetry oF matrices [a] and [b], Alsoy these are 
banded matrices, enabling one to store elements within 
one half the band width alone entailing thereby savings 
in space and computational times. 

Use was made of the NAG (Numerical Algorithms 
Group) FORTRAN library [?6] routine F04ACF for the solu- 
tion of the system of equations (2.49). The triangular 
ring element mesh is shown in Fig. 2,7. 



0 ’^ 


00 „ mf it© 












a-iAPTER 3 


3.1 RESULTS AND DISCUSSIONS 

The diffusion model and its finite element 
formulation that have been developed in the previous 
chapter are used to find solutions for global and micro- 
temperature Fields in the half space X > 0, subject to 
a temperature-pulse boundary condition. System sensiti- 
vities to non-linearities arising from temperature depen- 
dent thermal properties have been accounted for. Two 
sets of materials/ Graphite/ Epoxy and carbon/carbon com- 
posites have been considered as specific examples. A 
PEM computer pacbaoe is developed for this purpose and 
executed on the DEC 1090 svstem. This is attached as 
Appendix 1. The results have been presented in the 
graphical form and are discussed below. The temperature 
de'pendonce of thermal properties have been shown xn 
Appendices 2 and 3. 

The average values of the nodal temperatures at 
each cross-section of the fibre and the matrix have been 
plotted against the axial length of the Graph Ite/Epoxy 
composite. In Pig. 3.1. This has been done for a 
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particular Instant of bime just before the withdrawal of 
the constant temperature source at the boundarv. 
represents the fiber temperature and the matrix tem- 

perature. It is seen that the fiber is at a higher tem- 
perature than the matrix through a major length of the 
composite. It tends to the matrix or composite tempera- 
ture at distances away from the source. It can also be 
observed that the non-linearity of temperature properties 
provide lower temperature values in comparison to the 
linear case. This shows that the decrease in thermal 
conductivity with increase in temperature for Graphite 
plays a dominant role in the resulting temperature dis- 
tribution, 

' The average temperatures as explained above have 
also been plotted (Fig. 3.2) for the situation just after 
the withdrawal of the constant temnerature source in order 
to study the decay of the temperature fields in the two 
materials of the composite. Since the thermal conduc- 
tivity of tho matrix material is less than that of the 
fiber material, the temperature drop in the matrix mate- 
rial is much slower and consequently it is at higher 
temperature near the boundary. This also accounts for 
greater temperatures of the fiber at intermediate regions. 
It can be noticed that non-linearities provide for faster 
growth and decay of the temperature fields. 
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bohaviour of average temperatures and 

T 2 at an intermediate cross-section with the passage of 

shown in Pig. 3,3. evident, temperature builds 
up and docavs faster in l*e fiber than in the matrix 
accounting for the higher thermal conductivity of the 
fiber material. The figure also shows the difference 
b,^twoon the linear and non-linear behaviour of the com- 
po,‘; 1 to. 

Tho ability of tho analy.sis to prebict the tem- 
pera turo rnlcrostructure is illustrated in Fig. 3,4, The 
figure shows tho temperature fields across an intermediate 
cross-section a short while before and after the with- 
drawal of the constant temperature source. 

A short while before removal of the source, the 
temperature across the fiber is at a higher value than 
that across the matrix and is oractxcelly uniform. In 
the matrix core, it gradually reduces from a value of the 
fiber temperature at the interface to a minimum at the 
unit cell boundary. This is due to the lower thermal 
conductivity of the matrix material. 

A little while after removal of the source, 
however, the trend reverses. The fiber, having greater 
thermal conductivity, loses its heat faster and is conse- 
quently at a lower temperature than the matrix. 
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The material properties used in deriving the 
results discussed above represent a Graph it e/Fpo^y fiber 
reinforced composite. In this case, the thermal proper- 
ties depend rather weakly on temperature. The results 
however, imply that temperature dependence of thermal 
properties can significantly influence the results. In 
an effort to demonstrate this point, a similar initial 
value problem was posed and solved for a carbon/carbon 
composite. The material properties used in the evalua- 
tion arc provided in appendix III. The following dis- 
cussions reflect the observations made. 

Figure 3.5 represents the temperature distribu- 
tion along the axial length of the composite before with- 
drawal of the source. The results indicate that non- 
linear it, ies have significantly altered the distributions. 
It is worthwhile to note that temperatures calculated by 
taking into consideration the non-linear material beha- 
viour, which is representative of actual working condi- 
tions to a greater extent, are much lower than those 
obtained by the assumption of linear behaviour. A possi- 
ble consequence of this knowledge would perhaps lead us 
to employ lesser insulating thicknesses than we possibly 
would have, entailing cost savings. 
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On wi-t.hdrnwal of the constant temperature source 
the tompornture distributions are greatly altered as can 
be seen from f.he figure 3.6. The figure also provides 
ample evidence of the effect of non-linearities on the 
temperature distribution. It can be observed that 'hon- 
1 inear" temperatures are much lower than the " linear” 
temperatures for most of the length of the composite. 

Figure 3.7 shows the temperature time history at 
an intermedlnte cross-section of the composite. As long 
as tVie .<-!Ourc:e exists, the temperatures build up. On 
removal of the source however, there is a decay in the 
temperntxxr'o fields. The figure also illusbratos signifi- 
cant deviations caxised due to non-linearities in the 
thermal behaviour of the composite. 

The temperature distribution across an inter- 
mediate section, just before, and after withdrawal of the 
constant tc'mTH''rat.ure source are shown in figure 3.8. As 
expc'crtod, it is observed that the temperature builds up 
and df.'cavs faster tacross the cross— secticn of the fiber 
due to the higher thermal conductivity of the fiber mate- 
rial as compared to that of the matrix material. 
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3.2 CONCLUSIONS 

The beet conduction in uni~directional fiber 
reinforced composites subjected to a temperature pulse 
has been discussed in the present work. The effect of 
non-linearities on the temperature distributions have 
been amply demonstrated. The following conclusions may 
be drawn: 

(i) Non-linearities arising from temperature-dependent 
thermal properties significantly influence the 
temperature distributions. 

(ii) The non-linearities in case of both Graphite/ 

Epoxy and Carbon/Carbon, lead to temperature 
fields having magnitudes lower than those in 
the linear case.-. 

(iii) The non-linear analysis of the problem leads to 
an important result that less composite material 
is actually needed for insulation purposes than 
what is predicted by the linear solutions - 

(iv) Since the non-linear behaviour of the composite 
structures is a true representation of the 
materials used in practice/ it is recommended 
that non-linearities must be accounted for in 
the heat diffusion analysis. 
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APPENDIX IJ PROGRAM GISTTNS. 
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THIS IS A PROGRAM TO CALCULATE THE TEMPERATURE DISTRIBUTIOH 
IN COMPOSITE MATERIALS SUBJECTED TO HIGH TEMPERATURE PDLSES, 

THE PROGRAM USES THE FINITE ELEMENT METHOD OF SOLUTION » 

DIVIDING THE DOMAIN INTO A NUMBER OF FINITE ELEMENTS. 

THE PROGRAM USES THE r04ACF ROUTINE OF THE NAG [NUMERICAL ALGORl' 
THMS GROUP] LIBRARY TO SOLVE THE SIMULTANEOUS EQUATIONS BY 
CHOLESKY'S DECOMPOSITION TECHNIQUE. THIS ROUTINE SOLVES 
EQUATIONS OF THE TYPE AX«B, WHERE A IS A SYMMETRIC , BANDED 
MATRIX. ONLY THE LOWER TRIANGULAR ELEMENTS OF THE ^ 

COEFFICIENT MATRICES ARE USED, TO SAVE STORAGE SPACE AND SPEED 

» 

PARAMETER“ipIa8o7iP2=l2o7iP3=5liP4=7i7iBl=lO 

COMMON NN(IP2,3) ,XC0RD(IP1) ,Yc6rD( 1P1),MATL(IP2) ,DELTA,T0C IPl 3 
ItlTCIPi) 

DIMENSION AE(3,3) ,BEC3j,3),A(IP4,IB1),X(XP1) 

OPEN(UNITa50 rDEVICEs'DSK*’ 5 



GiNiRAL"it5FORHATi5N’" ' 

READfSO.lOOO) TITLE 

writecsLioodtitle 

read ( 50, ♦)NN0DE,NELH,NMAT, DELTA, ISPN, ITER, TMP 
WRITE CSi, 10023 NN0DE,NELM,NMAT, DELTA, ISPN, ITER, TMP 

DIMENsioNING’ERROR'^DETECTORS 
if(nnode.gt.ipi3Go to 5 
lF[NELM.GT,iP2)G0 TO 6 
IFCI5PN,GT,IP3)G0 TO 7 
GO TO 8 

WRITE{5,U13?STOP 

WRITE(5,n25;STOP 

WRITE(5yU3);STQP 

CONTINUE 


NODAL INFORMATION 
WR1TEC51.1003) 
NaNNODB^ISPN 
DO 10 K»1,NN0DE 


mm 

3$« 

j400 


uu i v , ivnuvci 

READ ( 50, ♦ 3 NDNO . YCORD ( K 3 I XCORD (K 3 , TO (K ) 
WRITE(;5l,lO04)K,XCORDCK),yCORDCK),T0(K) 

CONTINUE 

element'infoRmation"""’""'”*" 

WRITE(51,10053 
do 20 Jsl.NELM 

READ C50,#5eLN0,CNNI!JiI 3,1*1, 3) #MATL(J)^ 
WRITE(;5i,10063J,CNNca,l},I«lr3)»MATL(J) 

CONTINUE 

CALCULATE*’THE"BAND"wiDTH 

IBaO . ; 

DO 30 IEa4,NELM 

jBMAXO(NNaErl),NN(lE,23rNN(IE,33)**HINO(NN(IE,l),NNCIE/2) , 
1NNCIE43)3+1 

IF C IBi ,GT, IB ) GO TO 40 
WRITEC51,1145 IB ; STOP^ ^ ^ 

CONTINUE _ ^ 

GENiRATi''THE’’B7c7"cODE*FOR*THE NODES 

DO 50 Isl,NNODE 

11(13*0 

CONTINUE 

DO 60 Isl,ISPN 

READ(60,»5k 

ITCKIal : :■ . , , 

. eoNTiNUE: ■/ 

NOw’'RENUMBER"iHB"EQUAllONS’'RBMOVjNG^ONis*'i5lTH-CONST^ 

: K2»NNQDE+l ■ 

DO 80 Is,l,NNDDE 
lF(ITCI)*iQ.0)G0 T070 
K-2*K2M.l^ ■ : 

GO 'TO. :S0 :r 
Ki«Ki+i' 

iflllsKi ' : ; ■■ ■ 
eoMpuE. 

DO 90 lalvNNODE 
TtIT(I)3*TO(I) 
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CONTINUE 

DO 95 I=1,NN0DE 

TOdJaTCX) 


500 

C 

C 

111 

tl 2 

113 

111 

1000 

1001 

1002 


1003 

1004 

1005 

1006 
1007 
1006 


TIMEsO.O 

IRsl ; MaIB-1 t MisiB 


IjOdp for increments in time 
DO 500 NIT=1/ITER 
IFCTIHE.liT.TMP) GO TO 106 
DO 105 KaN+l,NNODE 


T(K)aO.O 

CONTINUE 


CONTINUE 
CONTINUE 
DO 100 1=1, N 
DO 100 J=1,IB 
Ad.JiaO.O 
CONTINUE 
DO 110 1=1, N 
TO(I)=Td) 

Td)s0.0 

CONTINUE 

DO 300 IE=1,NELM 
CALti CONDCIE.AE) 

CALL HTCAP(2E,BE) 

DO 200 Ilal,3 
IGalTCNNCIEdi)) ^ 

IFdG,GT,N)GO TO 200 
DO 200 Ji=l,3 
JGaITCNN(IE,Jl)) 

IFCJG.GT.N) GO TO 150 
JG1=JG-I,G + IB 

IFC JGl ,LE. 0 .OR, JGI .GT, IB ) GO TO 145 
ACIG,JGl5=AdG,JGn+AEdl,Jl)+BEdl,Jl) 

CONTINUE 

TdG)aTdG) + (BE(Il,Jl)-AE(Il,Jl))^TO(JG) 

GO TO 200 

TdG)aTdG)-2,*AEdl,01)»TOCJG) 

CONTINUE 

CONTINUE 

TIHEsTIME+DELTA 

CALL^F04ACF(A.IP4,T,IP1,N,M,IR,T,IP1,A,IP4,M1,IFAIL) 

IFCIFAIti.EQ.OlGO TO 555 

WRITEC5»*)IFAIL 

STOP 

CONTINUE 

WRITE(51 ,1007)TIME 

WRIXEC5lJlO08jdCIT(I)),Ial,NNODE) 

FORMAT^STATEMENTS"’”""" " 

FORHATCIH, 'INCREASE IPl!) 


FORHAX( 20 A 4 ) 

FORMATclH S0A4J 

FORMAXaH j^NUMBER OF N 0 DES='’, 13 .// 1 H ,'NUMBER OF ELEMKNTSf',, 

113. //iH .Number of materials =:. 12 . //ih .'incremental ,time 5 % 

2Fli,8,//iH /'NO, OF NODES WITH SPECIFIED TEMPS.a' , 13 .//IH , 'NO, 
lOP iTfeRATIo/fS-'Iia./ZlH ,'IHPULSE DURATlONa '.Fl2;8/^) 

15 X,' 0 '. 5 Xl'K'; 3 X.'MATpIAL'/y^^^^^ ^ - 

P0rAaT(6X,I3,3I6/4X,I4) ^ 

F0RMATC//1H ,'TiAe»%F11.5) 

FORHAT(S(lPEi6.43J 


V/IH ,7X>'ELE'y3Xi'I' 


CLOSECUNITsfSO) 
CLOSE (UN1T»51) 
END : 


6700 


SUBROUTINE CONDTN.AB) 

PARAMETER IP 1 ago 4P2»120„.^ 

COMMON NNCIP2> 3) rXC0RDdPl) ,yC0RD(IP1 
; 1, IT (IP d. ^ ; 

DIMENSION AE(3/3) 

FIND THE NODE AUMBERS OF THE ELEMENT 
'IaMNCN, 13 r:, 


aaaaaasSrasstssasaassaa&saisasaaaBaaaaaasaaaaBB'asBaasaBsiteaaseaiaseisa: 

SUBROUTINE TO CALCULATE ELEMENTAL CONDUCTIVITY MATRIX, 
SUBROUTINE cond;n,ab 3;.^ r 

COMMoi’’NN C ipivif rXCORDCl^ r YCORD ( IPl ) , M ATL ( IP 2 ) ^ DELTa ,T 0 (ipi 


K=NN(N,3) 

PIMf) THE COORDTNATES OP THE NODES 
XlsXCDPDCI) 

X2aXC0PI)CJ> 

X3=XCaPDfKT 

YlaYCORD(I) 

Y2-YCOPO(0) 

Y3=yCORDCK3 

FIND THE ELEMENT AREA 

(^l*^2-Y2*Xl) + (Y2*X3-X2*y3) + (Y3*Xl-Yl*X3n/2.0 

BlBy2"Y3 

B2aY3-Yl 

B3BY1-Y2 

C1«X3-X2. 

C2«X1-X3 

03»X2"*X1 

FIND .’THE CONpyCTIVITY AT THE^AVERAGE TEMP OF THE EfaEMENT. 


TAVs(TO(ITCI))+TO(ITCJ))+TO(1T(K)>)73;0 
CONDUC" “ 


SPECTFY THE CONDUCTIVITY-TEHPERATORE RELATIONS. 
TFCHATL(N).E0.2) go to 61 
C0HX=1069,42 
C0MYall0,84 

COMXs 1D69,42-1UR.78*CTAV/1000,)+646.40^((TAV/1000,)»»2J 
l^l79.327*aTAV/l500.)**3j+l9,i263*((TAV/1000,)**43 
CONYsllO,84-106,12*(TAV/lOOO,5+54.24»( fTAV/iOOO. )#»2) 
1-13.5?4((TAV/1000,)**3)+1,1647*( (TAV/1000,)**4) 

GO TO 71 
CONTINUE 
C0NX=178.2O 

C0NX«178, 20-31. 678*CTAV/1000.)+6,015»CCTAV/1000 , )**21 

CONYsCONX 

CONTINUE 

FOP CONVENIENCE DEFINE QUANTITIES XAVjDX.DY 

XAVaai+X2+X3)/3,0 

DXaC0NX>l<XAV/C4,0*AREA) 

DY3CONY*XAV/(4.0tAREA) 

ASa , 1 )sBl*DX*Bl+Cl*DY*Cl 
AE(l,2)sRl*DX*B2+Cl»nY*C2 
AEC1#3 )s 81»DX*B3+C1»DY»C3 
AEC2,2)aB2*DX>l'B2+C2*DY»C2 
AEC2,3)aB2#DX^B3+C2»DY»C3 
AEC3,35aB3*DX*B3+C3*OY*C3 
AE(2 l)rAE(l,2) 

AChyl JsAEd.aj 
AEU;2)aAEC2,3) 

RETURN 

END 

SSSSSaS!SSBBSBSB9SBBBSSSSBB=S«S3sasSSSSSBSsSs:saSSBBSSSaSt»SSasS 

SUBROUTINE TO CALCULATE THE ELEMENTAL HEAT CAPACITIES,: 
SUBROUTINE HTCAPCN.BE) ^ ^ ^ ^ 

PARAMETER IPlaBO , IP2«i20, lP3a5, IP4a75. IBlslO ^ v 

COMMON NNCIP2,3),XCORD(IPl)rYCORD(IPl5 ,MATL(IP2) ^DELTAiTOIlPl) 
l.ITeiPl) 

DIMENSION BE(3/3) 

GET THE NODE NUMBERS OF THE ELEMENT 

TaNN(N,n 

UaNN(N,2) 

KaNM(H,3) 

GET THF. CO-ORDINATES OF THE NODES 
XlaXCORDCI) 

X2aXCORD(J5 
X3bXC0RD(K) 

YlaYCORDCI) 

Y2ayC0RDCJ) 

Y3aYC0RD(K) ■ 

ARFAaABSC CY1*X2-Y2*X1) + 

ELUae ,0#Xt+2.0*X2+2.0* 

EL12s2,0*Xl+2.0»X2+X3 
EL13a2.0»Xl+X2+2.0*X3 
EL22a2,0»Xl+6,0*X2+2.0*X3 
EL23aXl+2,0*X2+2,0*xV^;, 

EL33a2,0*Xl+2.04«X2+6.0*X3 
FIND THE AVERAGE TEMP ^ . ,, 

TAVoCTOC ITd ) ) +Tnc IT f J) 3+T0( ITfK) ) )/3 
DEFINE THE HEAT CAPACITIES HERE 
IFCMATLCN).EQ.2)G0 to 81 


-Y31'X2) + CY3*Xi-Yl*X3))/2,0 


HCsO. 2937 
HCaC 


. rtA§;f 

GO 'Tb^ 91 ^ : 

CONTINUE V 

HcSov277l+l. 036 I^CTAV/IOOO ;)tO. 7984* t(TAV/ipOO^)»*3) +0.2^9^^ 



1TAV/I000.5**3)-0.0412*((TAV/I000,)>t!»4) 

CONTINUE 

C0ErF=HC*AREA/60, 

FACTb2, /DELTA 

RE( 1 , 1 5aCOEFF*RLl 1»FACT 

DECl,255C0EFF*RL12*rACT 

WECi ,3)5CQEFF*EL13*FACT 

BE(2,25sCOEFF*EL22*FACT 

RE(2,3lsCOEFF*EL23*FACT 

BEC3,3)3C0EFF*EL33*FACT 

REf2,n3BE(l,2) 

RE(3#n«BEa,3) 

BEC3,2)sBE(2,3) 

RETURN 



APPiaSIDIX 2 


Material properties for Graph it e/Epoxy composite 



Axial 

j,(a) 

XX 

Thermal Conductivity 

Heat capacity 

^«x) 


^Radial 

K<"> 

rr 

Thermal Conductivity 

Fiber 

(a=l) 

k'D = 

XX 

(l) 

10-0.01 T^^^ 

^^=0,875 ■K).004T^ 
-5xlO“®(T^ 

k(J> = 

rr 

{ t) 

8-O.OOS 

Matrix 

(a=2) 

j^(2) _ 

XX 

(2 ) 

1 + 0.005 T^ ^ 

4^^^=140.002 T^^^ 

+1.5xlO“^(T^^^ 

k< 2> = 
rr 

(2 ) 

1 + 0.005 T^*^ ’ 




APPENDIX 3 


Material properties of carbon/carbon composite 


f ^ « E ^1000^ 

n=0 



Bq 

'^1 

^2 

. ^3 

, ^4 

K<1>' 

XX 

1069.42 

- 1119.78 

646.40 

179.327 

19.1263 

^(2) 

^xx 

178.20 

31.67B 

6.015 

0 

0 

^rr 

110.84 

- 106.12 

54.24 

13.52 

1.1647 

k(2> 

rr 

173.20 

31.678 

6.015 

0 

o 


0.29 37 

1.0962 

- 0.3448 

0.3129 

-0.0436 


0.2776 

1.0 361 

- 0.7934 ■ 

0.2953 

-0.0 412 

Note : 


[k^"'] = 

L rr ^ 

ca 1 /cm~h r- ” C 




= cal/cm^-®C 
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